Single-cell quantification and dose-response of cytosolic siRNA delivery

Endosomal escape and subsequent cytosolic delivery of small interfering RNA (siRNA) therapeutics is believed to be highly inefficient. Since it has not been possible to quantify cytosolic amounts of delivered siRNA at therapeutic doses, determining delivery bottlenecks and total efficiency has been difficult. Here, we present a confocal microscopy-based method to quantify cytosolic delivery of fluorescently labeled siRNA during lipid-mediated delivery. This method enables detection and quantification of sub-nanomolar cytosolic siRNA release amounts from individual release events with measures of quantitation confidence for each event. Single-cell kinetics of siRNA-mediated knockdown in cells expressing destabilized eGFP unveiled a dose-response relationship with respect to knockdown induction, depth and duration in the range from several hundred to thousands of cytosolic siRNA molecules. Accurate quantification of cytosolic siRNA, and the establishment of the intracellular dose-response relationships, will aid the development and characterization of novel delivery strategies for nucleic acid therapeutics.

Multiple positions of HeLa cells stably expressing YFP-galectin-9 were imaged before treatment with lipoplexed AF647-siRNA (siRNA-to-lipid ratio 2:4 pmol:µL) followed by image acquisition every 5 min with an Airyscan confocal detector of half of the positions (illuminated). At the end of acquisition i.e. t = 6 h, non-illuminated cells were imaged again. Each cell was then manually inspected for recruitment of YFP-galectin-9 on vesicles containing AF647-siRNA lipoplexes. (a) Illustrative images of illuminated and non-illuminated cells imaged before (0 h) and after (6 h An algorithm for automated detection of siRNA release was used, to identify cytosolic release events above one of two prespecified rolling thresholds (see Methods for details). Representative cytosolic siRNA measurements from one cell is shown, with a detected release event at frame 37. Source data are presented in the Source Data file. Fig. 7 | Model estimation of siRNA release quantity using Bayesian information criteria. Three functions (exponential, constant and double exponential) were used to fit a curve to the measured siRNA concentrations during endosomal release, as described in Supplementary Note 1. A relatively noisy example is shown to illustrate the selection of the trajectory with the best fit using Bayesian information criteria (BIC). Model 2 (center panel) is selected given the lower BIC despite model 3 (right panel) having a higher R 2 (reflecting over-fitting). For each model, ten fitted trajectories using different (random) starting points for the function is shown (lines). Circles represent siRNA measurements. Source data are presented in the Source Data file. HeLa cells expressing d1-eGFP were imaged every 5 min during treatment with 50 mg mL-1 cycloheximide (CHX) or DMSO (control). The d1-eGFP expression is shown relative to t = 0 (start of acquisition immediately after addition of CHX or DMSO). D1-eGFP half-life (t 1/2 ) was estimated to ~48 min. Mean and 95% confidence interval (shaded area) is shown. N = 218 and 102 cells for CHX and DMSO, respectively, from three independent experiments. Source data are presented in the Source Data file. To control the number of release events, the volume of prepared siRNA-lipoplexes added to cells were adjusted to 5 µL, 10 µL or 16.7 µL. A confocal microscope with Airyscan detector was used for live-cell imaging, followed by single-cell analysis. The siRNA release magnitude in individual cells is shown for each condition. Box-plot shows centre at median, box bounds 25th and 75th percentiles, whiskers min and max. N = 166, 79 and 105 cells from five, one, or two independent experiments, respectively. Source data are presented in the Source Data file.    HeLa cells expressing d1-eGFP were imaged every 5 min with a confocal microscope. D1-eGFP fluorescence intensity was measured using masks based on cell boundaries (entire cells) or cell nuclei only, and is shown relative to t = 0 (start of acquisition). Line is mean, shaded are is 95% confidence interval. N = 131 cells from three independent experiments. Source data are presented in the Source Data file.

Supplementary Fig. 19 | Flow cytometery gating strategy for eGFP knockdown experiments.
Gating strategy used to select intact cells in d1-eGFP knockdown experiments, presented in Fig. 3a. Control cells were analysed first and a gate was created that served as a template for siRNA treated samples. N Here, we descri�e how the pea cytosolic si �� uorescence intensity is estimated �y modeling release events. �f o�servations were exact and continuous, the highest measured value of the si �� �uantity would have �een the most accurate estimate of the pea si �� concentration. However, since si �� is o�served with measurement error and only at discrete time points, we cannot o�tain this data directly from o�servations. �nstead we �uild a mathematical model for the data, which in turn provides an estimated maximal value of the si �� amount, which we denote siR A The data on which the model is �uilt is the following: For each cell, we have o�servations of the si �� uorescence intensity, at frames −� to � relative to the release event occurring at frame . Three different types of distinct trajectories of the change in si �� uorescence intensity are fitted to this data. The trajectories were selected �ased on the o�served inetics of si �� release and redistri�ution. First is the constant trajectory, where a sudden signal increase is o�served at the release, which then levels off to a constant �Fig. �f�. �econd is the exponentially decaying trajectory �Fig. �e�. Third is the case with two su�se�uent release events �dou�le exponential� �Fig. �g�.
hen examining si �� measurements manually, it is often possi�le to see which type of event has occurred. To have the model autonomously decide what trajectory is most appropriate, all three trajectories are fitted for each cell. The model with the �est fit is then selected using �ayesian information criteria ����� ���. The criteria will punish models using too many parameters, thus penali ing overfitting. This is re�uired as the dou�le trajectory will typically give a �etter fit than the single trajectory, ma ing simpler measures for evaluating goodness of fit �li e R 2 � disadvantageous, since it would always select the most complex model.
Here we mathematically define the three trajectories used. �ll three descri�e a fast linear growth after the first si �� release, then followed �y different �ehaviours. The trajectories are denoted f (t; θ) where the first argument is time and second argument is parameters.
The function which we fit is the following, consisting of two terms The second term descri�es the linear growth from the release event �at frame θ 3 ) until frame having the largest value, and the first term descri�es the trajectory after the growth is complete.

10
Here the three parameters represent θ 1 − time after which the linear growth is complete, θ 2 − the value after the linear growth is complete and θ 3 − time where the linear growth starts. For this trajectory siR NA = θ 2 .
The function which we fit is the following function consisting of two terms.
Here we have the same parameters as the constant function, with the addition of three more: θ 4 − logarithm of the largest value from which the exponential function decays, θ 5 − scale on which the exponential function decays, θ 6 − shape of the exponential function. For this function This function is identical to the exponential trajectory except that it has one additional growth with exponential decay after the initial event, at a time point which needs to �e estimated �y the model. Thus, in total it has ten parameters.
For each cell, functions are fitted using least s�uares ���� estimate: for the three trajectories f 1 (t; θ), f 2 (t; θ) and f 3 (t; θ). �s these functions are non convex �i.e. with potentially multiple local maxima�, ten different random starting points for θ is used for each function, and the parameter with the smallest �� value is selected. �xamples of the different types of functions found are illustrated in �upplementary Fig. ��. �fter finding the function with the �est o�tained fit for each trajectory type, we use ��� to select the �est of the three.
For the �est model, here denoted , we also calculate what proportion of the variance that is explained �y the fitted function, i.e. R 2 . This statistic is calculated as follows shown. �n this instance, the exponential trajectory had the lowest ��� and hence was chosen. �otice that the dou�le exponential trajectory had the �est fit �R �. However, from examination of the data there is no clear evidence of two release events. This is correctly identified �y the ��� and the simpler trajectory is selected. 7 7 �n order to examine the effect of the cytosolic si �� �uantity on e F noc down in single cells, an explicit function for the relation is needed. The function we used is the following: ds .
Here, t denotes time after the release event �occurring at time � and R is the si �� uorescence intensity. The reason for using an exponential is the function is non negative of the e F . The actual function inside the exponential is easier to understand in terms of its derivative, Here, the derivative ta es its minimum at time point �after the time of release� where the value is θ 1 R θ 2 �this function was found after numerical fitting�. �t the time of the release �time is ero� the value of the function is ero and also as time goes to infinity. The parameter θ 4 controls the scale on which the derivative changes.
Supplementary Note 2. Uncertainty analysis of siRNA release estimation We estimate cytosolic siRNA release amounts ( ) by fitting a function to reference concentration curve converted fluorescence measurements ( ) according to: Converted siRNA concentration measurements ( ) are related to uncorrected raw fluorescence measurements ( ) as, where, conversion factor from fluorescence to siRNA concentration (derived from reference measurements before each experiment) bleaching correction of releasing siRNA particle due to light exposure during imaging local background correction or zero-offset Uncertainty of the final siRNA concentration estimation ( ) will be dependent on the errors introduced in the estimations of each unknown parameter ( , , , ). The estimated relative uncertainties of these parameters can be expressed as: Relative error of the conversion factor. Estimated as the relative standard deviation of the reference measurements before each experiment. In average 7.8% (21 experiments).
Relative error of bleaching induced uncertainty. Estimated as the standard deviation of the number of frames each releasing particle is imaged before release, multiplied with the relative bleaching per frame. Values are derived using galectin-9 expressing cells (Fig. 3a, 3 experiments). 5.8 frames 0.7% 4.06%.
Relative error of the local background correction or zero-offset. Estimated as the standard deviation of the (mean) background intensity in the 10 frames before siRNA release (which is the window used for zero-offset correction). This value is divided by the final siRNA release estimation, to give the relative error on a per individual event basis.
Relative error of the model fit. Estimated by the least squares method to give the relative error on a per individual event basis.
Each error term is independent, where and are scaling factors and is additive in the model fitting. Thus, using the delta method the relative standard deviation of the estimated siRNA release is, The estimated relative errors of the individual release quantifications of siGFP-1 and siGFP-2 are presented in Supplementary Figure 12.